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ABSTRACT 

Accurate numerical solutions of the equations of hydrodynamics play an ever more im¬ 
portant role in many fields of astrophysics. In this work, we reinvestigate the accuracy 
of the moving-mesh code Arepo and show how its convergence order can be improved 
for general problems. In particular, we clarify that for certain problems Arepo only 
reaches first-order convergence for its original formulation. This can be rectified by 
simple modifications we propose to the time integration scheme and the spatial gradi¬ 
ent estimates of the code, both improving the accuracy of the code. We demonstrate 
that the new implementation is indeed second-order accurate under the L 1 norm, and 
in particular substantially improves conservation of angular momentum. Interestingly, 
whereas these improvements can significantly change the results of smooth test prob¬ 
lems, we also find that cosmological simulations of galaxy formation are unaffected, 
demonstrating that the numerical errors eliminated by the new formulation do not 
impact these simulations. In contrast, simulations of binary stars followed over a large 
number of orbital times are strongly affected, as here it is particularly crucial to avoid 
a long-term build up of errors in angular momentum conservation. 

Key words: methods: numerical, hydrodynamics, galaxy: formation 


1 INTRODUCTION 

Gravity and hydrodynamics provide the foundations for de¬ 
scribing almost all phenomena in astrophysics. Often, the 
dynamics and geometry are sufficiently complicated that 
analytic solutions, especially for hydrodynamics, cannot be 
obtained. These limitations can be overcome by numerical 
approaches, which have developed over the recent decades 
into powerful tools for studying complex hydrodynamical 
flow problems. It is however an ongoing and persistent chal¬ 
lenge to derive ever better discretization schemes that re¬ 
duce numerical discretization errors to a minimum, while 
at the same time offering high adaptivity to different time- 
and length-scales and a high degree of parallel scalability. 
Astrophysical code development is therefore best viewed as 
an iterative effort which never is fully complete and fin- 
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ished. Rather, new generations of numerical ‘instruments’ 
(i.e. codes) should ideally improve their accuracy and speed, 
somewhat similar in spirit to the advances regularly realized 
in observational instrumentation. 

Traditionally, stationary mesh codes with or without 
adaptive mesh refinement (e.g. Fryxell et al. 2000; Teyssier 
2002; Almgren et al. 2010; Stone et al. 2008; Bryan et al. 
2014) and smoothed particle hydrodynamics (SPH) codes 
(e.g. Wadsley et al. 2004; Springel 2005) have dominated 
astrophysical fluid dynamics. Moving meshes and hybrid 
techniques combining Eulerian and Lagrangian meshes have 
been applied only rarely in astrophysics and related fields 
(e.g. Hirt et al. 1974; Gnedin 1995; Pen 1998). However, in 
particular in astrophysics they have not seen widespread use 
so far. 

Recently, moving-mesli techniques based on a Voronoi 
mesh have been proposed (Springel 2010; Duffell & Mac- 
Fadyen 2011), which offer a quasi-Lagrangian description 
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that inherits some of the advantages of SPH while retain¬ 
ing the accuracy of a Eulerian mesh-based description. These 
new methods have matured into interesting alternative codes 
that have already been widely applied in cosmological simu¬ 
lations of structure formation (e.g. Vogelsberger et al. 2012; 
Sijacki et al. 2012; Marinacci et al. 2013; Vogelsberger et al. 
2014; Genel et al. 2014), as well as in first star formation 
(e.g. Greif et al. 2011, 2012; Smith et al. 2014), and to a 
smaller extent in stellar astrophysics (Pakmor et al. 2013) 
and in problems related to planet formation (Duffell & Mac- 
Fadyen 2012; Munoz et al. 2014, 2015). It has become clear 
that it can be very worthwhile to accept the technical com¬ 
plications introduced by a dynamic and unstructured mesh 
given the benefits realized by combining most of the advan¬ 
tages of fixed mesh codes (e.g. better convergence for smooth 
flows, good shock capturing) and SPH techniques (e.g. in¬ 
trinsic adaptivity, lack of advection errors, small numerical 
diffusion). 

A well-defined and rapid convergence rate is a partic¬ 
ularly important feature of mesh codes, which are usually 
constructed to be at least second-order accurate for smooth 
flows, i.e. in the absence of shocks or other discontinuities. 
For SPH codes, in contrast, it is very hard to demonstrate 
proper convergence in the first place (Zhu et al. 2015). In 
practice, the convergence rate dictates how rapidly a nu¬ 
merical solution improves as more computational effort is 
invested, an aspect that is arguably even more important 
than the absolute size of the error itself. For example, while 
SPH often obtains a qualitatively correct result with accept¬ 
able error at low resolution, its poor convergence rate does 
not allow it to reach comparable accuracy to mesh codes at 
comparable numerical cost once the resolution is high. In 
contrast state of the art mesh-based codes are expected to 
show at least second order convergence. 

Moving-meshes in general have been studied extensively 
(see, e.g. Thomas & Lombard 1979; Guillard & Farhat 2000). 
In this paper, we reexamine the accuracy and convergence 
rate of the Arepo code (Springel 2010), which is at present 
the most widely employed implementation of the moving- 
mesh technique in astrophysics. We address weaknesses in 
the original version of this code and show how they can 
be overcome to improve the overall accuracy of the scheme. 
This directly affects the range of applicability of the code 
and is hence important for further maturing the moving- 
mesh approach in general. 

The paper is structured as follows. In Section 2, we con¬ 
cisely review the essential parts of the implementation of 
moving-mesh hydrodynamics in Arepo and analyse its con¬ 
vergence. In Sections 3 and 4, we describe improvements to 
the time integration scheme and the gradient estimate, re¬ 
spectively. We apply the new implementation to test prob¬ 
lems in Section 5, show results for realistic applications in 
Section 6 , and discuss the implications of our work in Sec¬ 
tion 7. 


2 MOVING-MESH HYDRODYNAMICS IN 
AREPO 

The moving-mesh code Arepo solves the Euler equations us¬ 
ing the finite-volume approach on an unstructured Voronoi 
mesh that is generated from a set of points (Springel 2010). 


At any given time during the simulation, the mesh can be 
uniquely constructed given only the positions of the mesh¬ 
generating points. Since the mesh-generating points move, 
the geometry of the mesh and its connectivity change over 
time. 

The Euler equations can be written as a system of hy¬ 
perbolic conservation laws for conserved quantities U and a 
flux function F(U), 


9U 

dt 


+ V ■ F = 0. 


(1) 


For the standard Euler equations the conserved quantities 
are mass, momentum and energy. They define a vector of 
conserved quantities and an associated flux function as 
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where p, v, P, and e are density, velocity, pressure, and 
total specific energy of the fluid, respectively. The latter is 
defined as the sum of the specific internal energy u and the 
specific kinetic energy |v J , thus e = u + |v 2 . The system 
of equation is closed by an equation of state P = p(y — l)e 
with the adiabatic index 7 . 

To solve the Euler equations, the state of the fluid is 
discretised using the cells of the Voronoi mesh. To this end, 
averages of the conserved quantities U of the fluid are com¬ 
puted for each cell through integration over the finite volume 
of a cell i, 


Q, = / UdV. (3) 

JVi 

yielding a census of the conserved physical quantities in the 
cell. This state is then evolved in time by 

Q" +1 = Q" - At ^ AfjF” +1/2 , (4) 

i 

where Aij is the oriented area of the face between cells i and 
j and F is a time-averaged approximation of the true flux 
Fy across the interface between cells i and j. 

The calculation of the fluxes F is done in the original 
version of Arepo using a MUSCL-Hancock scheme (van 
Leer 1984; Toro 1999), which has been shown to provide 
second-order accuracy in time and space for all kinds of sta¬ 
tionary meshes (see, e.g. Fryxell et al. 2000). This method 
uses a slope-limited piece-wise linear spatial reconstruction 
step in each cell and a first order time extrapolation of the 
fluid states by half a timestep to obtain the states on both 
sides of all interfaces. Finally, a Riemann solver uses the 
states on both sides of an interface to compute the flux that 
is exchanged during the timestep. 

The extrapolation is carried out in primitive variables 



( 5 ) 


which are straightforwardly calculated from the conserved 
quantities and the geometry of a cell with the equation of 
state, combined with estimates for their local spatial gradi¬ 
ents. By construction, they represent the values at the center 
of mass of a cell. Then, the left and right states of an inter¬ 
face are computed by a linear spatial extrapolation from the 
center of mass s of a cell to the center f of a cell face, and 
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by a half-step prediction forward in time (where At is the 
full timestep) as 


„ r/ dW 

W L r = Wl,r H—x— 

dr 


L,R 


,, s , dW 

(f “ Sl ' r) + ~Bt 




The time derivatives of the primitive variables in 
Eqn. (6) are not calculated directly. Instead, the Euler equa¬ 
tions are used to express them in terms of the primitive 
variables and their spatial derivatives only: 


dp 

dt 

dv 

dt 

dP_ 

dt 


vVp — pVv, 

(7) 

VP nr 

-vvv , 

P 

(8) 

qPVv — vVP. 

(9) 



The extrapolation and computation of the fluxes over an 
interface are all done in the rest frame of the (moving) 
interface, i.e. the interface velocity is subtracted from the 
equations above. This allows solutions that are manifestly 
Galilean-invariant, an important conceptual advantage over 
traditional Eulerian methods where the numerical trunca¬ 
tion error grows with the fluid velocity. 

In the original implementation of Arepo, the spatial 
gradients are calculated using an improved Green-Gauss es¬ 
timator that makes use of certain mathematical properties 
of the Voronoi mesh. Specifically, the gradient estimate for 
a primitive variable <f> in cell i is given by 


W>.-K I> 


A+AlE2 

T a A Tij 


where 
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is the center of mass of the face between i and j, and 
Tij = n — r j is the difference between the positions of the 
mesh-generation points of cells i and j, and r,j = |r;j| is its 
length. The sum runs over all cells that share an interface 
with cell i. Note that this estimate of the spatial gradient of 
a primitive variable depends only on the values of the prim¬ 
itive variable in the cell itself and its direct neighbours, and 
on the local geometry of the cell. If the values of the primitive 
variables sample an underlying smooth held at the positions 
of the mesh-generating points, then the gradient estimate is 
second-order accurate for an arbitrary mesh geometry, i.e. it 
reproduces a constant gradient in the held to machine pre¬ 
cision. However, the primitive variables Q used to compute 
the gradient estimates represent volume averaged quantities 
rather than the value of the huid at the mesh generating 
point of a cell. 

In the original AREPO implementation of Springel 
(2010), all computations that require geometrical quantities, 
e.g. volumes, areas, or positions of face centroids, use the ge¬ 
ometry of the mesh at the beginning of the timestep. Never¬ 
theless, it has been argued that the scheme is second-order 
accurate, based on simple sound-wave tests (Springel 2010) 
using the L 1 norm, and based on the much more-demanding 
stationary isentropic vortex how (Yee et al. 2000), where the 
error decreases quadratically with increasing spatial resolu¬ 
tion (Springel 2011). 


Figure 1. L 1 norm of the density field as a function of the number 
of resolution elements per dimension at t = 10 for the isentropic 
Yee vortex, simulated with the standard Arepo implementation 
using MUSCL-Hancock time integration and the improved Green- 
Gauss gradient estimates. 


In this paper, we define the L p norm of a continuous 
field as 



( 12 ) 


where / (r) is, for example, the deviation of the density field 
from its analytical value at a position r. For our finite volume 
discretisation this becomes 

i/p 


L p = 
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The L 1 norm of the density field for the Yee vortex (for 
the detailed setup, see Appendix A) evolved until t = 10 
with the original AREPO code is shown in Figure 1. In 
contrast to previous claims, the error only decreases lin¬ 
early with resolution, i.e. AREPO shows first order conver¬ 
gence only. The difference to previous results (in particular 
Springel 2011) arises from an inappropriate definition of the 
L 2 -norm employed there. 

On static Cartesian meshes, we find that Arepo is fully 
second order convergent for smooth problems. On moving 
meshes, however, an inaccuracy is introduced since the time 
integration only uses the geometry of the mesh at the be¬ 
ginning of the timestep, ignoring changes of the mesh ge¬ 
ometry during a timestep. One possibility to account for 
these changes in the time integration is to do an additional 
mesh construction at the middle of the timestep and use 
the geometrical properties at t = to + At/2 to calculate 
the fluxes. However, this requires at least one further (ex¬ 
pensive) mesh construction per timestep, nearly doubling 
the computational cost. Such an approach combined with 
a Runge-Kutta time integration scheme is followed in Tess 
(Duffcll & MacFadyen 2011). 


3 RUNGE-KUTTA TIME INTEGRATION 

Alternatively, one can abandon the MUSCL-Hancock 
scheme and try to use a different Runge-Kutta time inte- 
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gration scheme that avoids a mid-step mesh construction. 
Of particular interest for us is Heun’s method, which is a 
second-order Runge-Kutta variant that calculates the flux 
as an average of the fluxes at the beginning and the end of 
the timestep. Applying Heun’s method to our moving mesh 
leads to the following update equations of the conservative 
variables for a cell i 


q ; = Qr-At^^F^u"), 


Q 


r 

n+1 


n+1 


(14) 

(15) 


n i a i n 

r + At w 

q: - 
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r + — (w -(- w 


(17) 


Here, the vectors r denote again the coordinates of the 
mesh-generating points, w their velocities, and F,j is an 
approximation for the fluxes of the conserved variables over 
the interface between cells i and j. In principle it is possi¬ 
ble to allow the velocities of the mesh-generating points to 
change during a timestep (Duffell & MacFadyen 2011), but 
we choose to always keep them constant over the course of 
a full timestep. 

Note that this update rule also requires the geometry 
of two different meshes (r n and r'), and the fluxes have to 
be computed twice for every timestep. However, since for 
Heun’s scheme and a constant velocity w' = w" of the mesh¬ 
generating points over the whole timestep we can show that 


n + l n . At /n. f\ n . \ . n / 

r = r -|—— (w + w ) = r + At w = r , 


(18) 


we can reuse the mesh we constructed for the second half of 
the current timestep for the first half of the next timestep. 
Thus, we only need to construct the mesh effectively once 
per timestep, which keeps the mesh construction effort equal 
to Arepo’s original implementation. The fluxes need to be 
calculated twice as often, but the computational cost re¬ 
quired for this is a relatively small effort compared to the 
3D mesh construction needed in a moving-mesh code, and 
therefore does not increase the overall computational cost 
significantly. 

In practice, computing the intermediate state Q' 
through Eqn. (14) is inconvenient for several reasons. In par¬ 
ticular, it significantly complicates the internal bookkeeping 
of the conserved variables since we need to know Q” and Q' 
at the same time in Eqn. (16). While this can be resolved in 
principle, further complications arise when one tries to con¬ 
sistently implement this update scheme for individual time 
steps that vary locally. 

These problems can be circumvented in a simple way by 
extrapolating to the intermediate step using spatial deriva¬ 
tives, like in the previous MUSCL-Hancock scheme. This 


leads to an update of the conservative variables through 


q; 


W' 

___ n . C/ VV 
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(19) 

( 20 ) 
( 21 ) 


( 22 ) 


Here, we only need to do the time extrapolation for the 
primitive variables, and we again use the Euler equations to 
replace time derivatives with spatial derivatives. Note that 
this update scheme can be easily generalized to local individ¬ 
ual time steps. For the calculation of the fluxes between two 
cells on different time steps, the time extrapolation (Eqn. 19) 
is in this case done for each cell always from the last time the 
cell was active. At this time the estimates for its primitive 
variables and gradients were obtained. Similarly, the linear 
spatial extrapolation to the current center of the interface 
is always done from the center of mass for which primitive 
variables and their gradients have been calculated. Note, 
that with the time-extrapolation to the intermediate step 
our scheme is not a Runge-Kutta scheme anymore, but be¬ 
comes a hybrid between Runge-Kutta and MUSCL-Hancock 
schemes. 


4 LEAST SQUARE GRADIENT ESTIMATES 

Another source of accuracy degradation in the original 
Arepo implementation lies in the quality of the gradient es¬ 
timates, which can suffer in general problems. The issue ap¬ 
pears when a cell is distorted and its mesh-generating point 
deviates significantly from the position of the center of mass 
of the cell. In this case, the Voronoi-optimized Green-Gauss 
estimate of Eqn. (10) introduces a systematic error, because 
it assumes that the primitive variables are known at the 
mesh-generating points, whereas, by construction, the primi¬ 
tive variables represent volume-averaged quantities that rep¬ 
resent the value of the underlying field at the cell’s center-of- 
mass. Although the small steering motions added by Arepo 
to the velocities of the mesh generating points usually man¬ 
age to keep the mesh nicely regular with only a small offset 
between the mesh-generating point and the center-of-mass 
of a cell, the typical distance between these two points can 
amount to a few percent of the radius of the cell, enough to 
significantly degrade the accuracy of the gradient estimate 
in some situations. 

We investigate this in Fig. 2 explicitly, where we show 
that such a small deviation already compromises the accu¬ 
racy of the gradient estimate, spoiling its convergence rate 
for higher resolutions. In fact, a one percent deviation of the 
mesh-generating points from the centers of mass in an al¬ 
most Cartesian mesh already stops the improvement of the 
gradient estimate with resolution for resolutions better than 
320 2 cells. Even worse, for a completely random mesh in 
which the mesh-generating points are a Poisson sample of 
the simulation domain the gradient estimate does not im¬ 
prove with resolution at all. Note, however, that the hydro 
scheme still converges with first order, even though the error 
in the gradient estimate is constant. 
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N 


Figure 2. L 1 error norm of the Voronoi-optimized Green-Gauss 
gradient estimate of the density field for different types of meshes 
for the initial state of the Yee vortex at t = 0. 



Figure 3. L 1 norm of the least-square-fit gradient estimate of 
the density field for different types of meshes for the initial state 
of the Yee vortex at t = 0. 

To more robustly reach second-order convergence of the 
hydrodynamical scheme, we therefore need better gradient 
estimates that are correct to at least first order for severely 
distorted meshes where center-of-mass values are known, 
equivalent to requiring that linear gradients are still repro¬ 
duced exactly in this case. 

Assuming that we know the local geometry of the mesh 
and the values of a quantity <j> at the centroids of the cells, 
we can obtain such a gradient estimate through a local 
least squares fit. We note that polynomial least squares re¬ 
constructions are well known for unstructured meshes, and 
form a standard technique in particular for high-order ENO 
and WENO schemes (e.g. Ollivier-Gooch 1997). Also, such 
methods are well known for the estimation of gradients (e.g. 
Maron & Howes 2003). They have also been recently used 
for SPH and in new variants of mesli-less methods (Hop¬ 
kins 2014) and to compute the cell-centered magnetic fields 
in constrained transport schemes for unstructured meshes 
(Mocz et al. 2014). The method assumes that the quantity 



Figure 4. L 1 norm of the density field for the evolved Yee vor¬ 
tex at t = 10 when different schemes for the hydrodynamics are 
employed. For the detailed setup of the Yee vortex problem see 
Appendix A. 

<f) can be approximated everywhere within cell i as 

<A(r) »0(si) + (V<A) i (r-s i ). (23) 

The value of (f> at the center of mass of the cell r/>( s i) = 
4>i is known as a prerequisite, and we are now looking for 
an estimate of the gradient (V0) i . We also know (j> at the 
centres of mass of all neighbouring cells j, so we can require 
that our gradient estimate reproduces those 4>j as well as 
possible. Thus, for every neighbouring cell, we would like to 
have 

</>j = + (V0)i (s j - Si). (24) 

For N neighbouring cells this is an overdetermined set of 
N equations with only up to three free variables (or two in 
2D). We select the best linear approximation for (V</>) i by 
requiring that it minimises the sum of the deviations for all 
neighbours, 

Stot = 9j (<fe - <t>i ~ {^<t>)i (s j - Si)) 2 . (25) 

3 

Here, gj is the relative weight for neighbor j. Different 
choices for these weights are possible, but according to our 
experiments the results are not very sensitive to the par¬ 
ticular choice made. We follow de Vasconcellos & da Silva 
(2014) and set it to gj = Aij/ |s 3 - — s;| 2 . 

To minimise Stot, we use the normal equations which 
yield 

gjKji ® n V I S 3 - Si| 2 = 

j 

(& - 4>i) |sj - Si I , (26) 

3 

where riji = (sj — Sj) / |sj — s;|. This corresponds to a 
2 x 2 or 3 x 3 matrix inversion problem in two- or three- 
dimensional problems, respectively, which can be solved by 
an appropriate solver to obtain (V0) i . 

The accuracy of the new least squares fit gradient es¬ 
timate is shown in Fig. 3. On the nearly Cartesian mesh 
examined earlier, the gradient estimate is as accurate as 
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the Voronoi-optimized Green-Gauss estimate. For increas¬ 
ing deviations from the Cartesian mesh, the accuracy of the 
least squares gradient estimate slowly deteriorates, but im¬ 
portantly, it never becomes worse than first order. In partic¬ 
ular, even for the random mesh it still accurately recovers the 
linear gradient, unlike the Voronoi-optimized Green-Gauss 
estimate. 


5 TEST PROBLEMS 


To examine the difference between the original Arepo im¬ 
plementation that employs a MUSCL-Hancock time integra¬ 
tion scheme (MH) and the Voronoi-adjusted Green-Gauss 
gradient estimate (GG) versus the Runge-Kutta time inte¬ 
grator (RK) combined with the least squares gradient es¬ 
timator (LSF), we compare these schemes for several rep¬ 
resentative test problems. In obtaining the results for the 
different implementations, we always use identical code con¬ 
figurations, parameters, and initial conditions, and only vary 
the time integration method and/or gradient estimate. 

All test problems in this paper are run using a Courant 
factor of 0.3 to determine the timestep and employ the same 
strategy and parameters to evolve the mesh while keeping it 
regular. To regularise the mesh, the mesh-generating points 
are moved towards the center of mass of their cell. Their 
total velocity is calculated as 


Vvertex — V 


1 . VP 

-At - h v reg 

2 p 


(27) 


where v, P, p, At, and v reg are the fluid velocity, pressure, 
and density, timestep of the cell, and a regularisation com¬ 
ponent that is added purely to improve the mesh. It is given 
by 


v reg 


= emax (c, r | V x v|) 


d 

W\’ 


(28) 


Here, c is the sound speed in the cell, r the approximate 
extent of the cell calculated from its area (volume) in 2D 
(3D) assuming the cell is a circle (sphere), and d = r — s is 
the separation between mesh-generating point and center of 
mass of the cell. Note that this separation is interesting in 
its own right, as it can be used to quantify the regularity of 
the mesh. Moreover, e is defined as 


e = max 


^0, min 


^0.5, 0.5 


a — 0.75/A \ 
(125/3 )) 


and a is given by 


a = max 

i 


( l 2 Aj \ 
V-^dim |r - r;| ) 


(29) 


(30) 


where the maximum runs over all neighbours of a cell and 
their interfaces connecting a cell to them. We choose the 
parameter fi to be 2.25 in all simulations in this paper. 



Figure 5. Conservation of total angular momentum for the Yee 
vortex at t = 10 when different hydrodynamical schemes are used. 
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5.1 Yee vortex 

We first reevaluate the convergence rate for the isentropic 
Yee vortex. Since this vortex flow is completely smooth, we 
in principle expect to see second-order convergence in the 
strong L 1 norm if the implementation is correct. As shown 
in Fig. 4, we indeed obtain second-order convergence up to 
very high resolution for Arepo in the moving-mesh case 


Figure 6. Density error compared to the analytical solution (top 
panel) and relative separation between mesh-generating point and 
center of mass (bottom panel) at t = 100 for the RKLSF Yee 
vortex. 
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if we use both the least squares gradient estimate and the 
Runge-Kutta time integration scheme. Using only one of the 
two improvements while keeping the MUSCL-Hancock inte¬ 
gration scheme or the improved Green-Gauss gradient esti¬ 
mate, respectively, does not change the convergence rate at 
all and makes the implementation drop back to first-order 
convergence in the L 1 norm. The detailed setup including 
the initial mesh can be found in Appendix A. 

It is interesting to note that the Tess code has been 
shown to be second-order accurate for an isentropic sound 
wave with a Runge-Kutta time integration scheme (Duffcll & 
MacFadyen 2011), while the same Voronoi-optimized Green- 
Gauss gradient estimate as implemented originally in Arepo 
has been used. With this gradient estimate, the convergence 
order is however expected to break down for meshes in which 
the mesh-generating points are not close to the centers of 
mass of their cells. The difference here is most likely that 
the one-dimensional propagation of the wave in this par¬ 
ticular problem only stretches the mesh, without displacing 
the mesh-generating points from the centres of mass of their 
cells. As discussed before, in this special case the Voronoi- 
optimized gradient estimate is sufficient and does not nega¬ 
tively impact the convergence properties. For the more de¬ 
manding Yee vortex this is not the case, however, and the 
convergence is expected to break down in Tess just as in 
the original Arepo code. 

Another interesting observation to make is that second- 
order convergence in Arepo for the new implementation as 
measured by computing the L 1 norm between the current 
density field and the analytical solution is only reached at 
very high resolution when the initial mesh-generating points 
are set up in rings around the center of the vortex. If an ini¬ 
tially Cartesian mesh is adopted instead, the convergence de¬ 
grades to first order at sufficiently high resolution. As shown 
in Fig. 6 there is only a weak correlation between the density 
error in a cell and its distortion, as both are largest close to 
the center of the vortex. 

This is fundamentally caused by the discretization of the 
analytical problem onto our mesh. There are two discretiza¬ 
tions involved, a first discreziation on the initial mesh to 
generate the discretized initial conditions and a second dis¬ 
cretization on the current mesh when we measure the L 1 
norm. If the structure of the mesh changes systematically 
between the initial and the current mesh, there is also a 
systematic difference in the two discretizations, which turns 
out to lead to a first order error that dominates the to¬ 
tal error at sufficiently high resolution. Because an initially 
Cartesian mesh will eventually be transformed to a spheri¬ 
cal mesh by the mesh regularization, the measured L 1 norm 
after the mesh changed will not be better than first order at 
high resolution. To overcome this problem, we arrange the 
mesh-generating points in the initial mesh already on circles 
around the center of the vortex. This mesh configuration is 
stable for the problem, thus there is no systematic change 
in the discretisation over time. 

Another approach is to look at error measures intrinsic 
to the discrete problem that do not require a second discreti¬ 
sation of the analytical problem. Analysing the conservation 
of total angular momentum (which also may be viewed as 
a norm) instead of the L 1 norm of the total density field is 
therefore considerably simpler and also more robust. In par¬ 
ticular, here the measured convergence rate is independent 


of the initial mesh, and the rearrangement from an initially 
Cartesian to a circular mesh does not give rise to additional 
errors. As shown in Fig. 5, the angular momentum conserva¬ 
tion is again only improved if both, a better time integration 
method and a more accurate gradient estimate, are used. In¬ 
terestingly, the angular momentum even shows third order 
convergence in the new implementation, although this may 
be the result of the special spherical symmetry of the vortex 
flow and may be lost for more general problems. 

5.2 Keplerian disc 

Evolving a cold Keplerian disc for many orbits is a com¬ 
mon problem in astrophysics, with applications from plane¬ 
tary to galactic discs. Using such a set-up in the limit of a 
pressure-less gas disc is a demanding test problem (Cullen 
& Dehnen 2010; Hopkins 2014). In fact, as highlighted by 
Hopkins (2014), many state-of-the-art static mesh or SPH 
codes have severe problems in coping accurately with this 
situation. Instead, the disc is typically destroyed during the 
first 10 orbits by these methods. 

The original AREPO code with the MUSCL-Hancock 
time integration and improved Green-Gauss gradients al¬ 
ready does quite well on this problem, as shown in Fig. 7 
where we display the disc after roughly 15 inner orbits with 
a resolution of 320x320 cells. However, whereas the outer 
parts of the disc are very stable, the inner boundary of the 
disc moves outwards with time owing to systematic errors 
in angular momentum conservation. Our new implementa¬ 
tion with Runge-Kutta time integration and least squares 
gradients works significantly better, keeping the disc essen¬ 
tially perfectly stable until t = 100, and showing similar 
errors only much later, at times t ~ 600 or later. These 
residual errors can now be much more efficiently improved 
to essentially arbitrary precision by an increase of the spatial 
resolution, thanks to the improved convergence order. This 
is in line with previous results on conservation of angular 
momentum in AREPO for more realistic problems (Munoz 
et al. 2015). 

These results compare very favourably not only to sta¬ 
tionary Eulerian grid codes on Cartesian meshes, but also 
to the new mesh-free hydrodynamical method proposed by 
Hopkins (2014). Note, however, that the optimal mesh con¬ 
figuration for this specific problem is most likely a polar 
grid. 


6 TESTING ON REALISTIC APPLICATIONS 

6.1 Evolution of a stellar binary 

In light of our previous results, it is interesting to investigate 
whether the improvements in the Arepo code proposed here 
affect the results of real world scientific applications of the 
code. Compared to the Yee vortex and the Keplerian disc, 
which are smooth hydrodynamics-only problems, or hydro- 
dynamical problems with a constant external gravitational 
field, one example of a more complex problem including 
discontinuities, shocks, and self-gravity is the inspiral and 
merger of a close binary system of two white dwarfs. 

Such simulations are essential to understanding the 
complex evolution of the merger process and to determining 
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Figure 7. Evolution of the surface density of a two-dimensional cold Keplerian disc. The top panels show the initial conditions and 
the evolved state at t = 100 obtained with the original Arepo implementation based on the improved Green-Gauss gradient estimate 
and the MUSCL-Hancock time integration. The bottom panels show the new implementation with least squares gradient estimate and 
Runge-Kutta time integration at f = 100 and t = 600, respectively. At the latter time, the inner disc has finished more than 250 orbits. 
The detailed setup can be found in Appendix B. 


the fate of the binary system (see, e.g. Zhu et al. 2013; Pak¬ 
mor et al. 2013; Dan et al. 2014). The initial mesh for the 
white dwarfs is constructed from shells that are tesselated 
using the HEALPIX algorithm to obtain regular cells with 
roughly the same mass (Pakmor et al. 2012). The simula¬ 
tion box is then filled by a low resolution Cartesian back¬ 
ground grid. Figure 8 shows the total angular momentum for 
the merger of two white dwarfs with masses of 0.65Mq and 
0.625Mq, an initial orbital period of 44s, identical to the 


setup of the same system in Zhu et al. (2013). We employ a 
mass resolution of 10 -6 Mq (high resolution) and 5 10~ 6 Mq 
( low resolution), respectively. 

The binary merges at around t = 200 s, after initial 
mass transfer reduced the separation and lead to runaway 
mass transfer. It then forms a differentially rotating merger 
remnant. The differential rotation in the merger remnant 
is crucial for its further evolution. As shown in Figure 8, 
conservation of angular momentum during the further evo- 


MNRAS 000, 1 11 (2015) 




Improving the convergence properties of the moving-mesh code AREPO 9 



f[s] 

Figure 8. Conservation of total angular momentum for the 
merger of two white dwarfs with an initial orbital period of 44s 
for the MUSCL-Hancock time integration and improved Green- 
Gauss gradient estimate and the Runge-Kutta time integration 
combined with the least squares gradient estimate. 

lution of the merger remnant is significantly violated for the 
implementation with MUSCL-Hancock time integration and 
improved Green-Gauss gradient estimates. The system loses 
a significant fraction of the total angular momentum present 
in the simulation. In contrast, the new implementation with 
Runge-Kutta time integration and least squares gradient es¬ 
timates conserves the total angular momentum in the sim¬ 
ulation to a relative error of about one percent even after 
many orbits for the high resolution simulation and still loses 
less than 10% of the total angular momentum in the low 
resolution simulation. 


6.2 Cosmological zoom simulation of galaxy 
formation 

An even more complex problem are simulations of galaxy 
formation and evolution. They not only involve self-gravity, 
but also feature high Mach numbers and turbulent flows, as 
well as a large number of additional source terms to model 
effects like radiative cooling of gas or the energy injection 
of evolving stars. In addition, sometimes explicit sub-grid 
models for unresolved physics are used that are introduced 
as a modified equation of state or pressure floors and the 
like. 

To compare the performance of our new code with the 
original Arepo version in this regime we have repeated one 
of the recent Milky Way galaxy formation simulations of 
Marinacci et al. (2013) and Pakmor et al. (2014). These 
are advanced calculations that can serve as an example for 
current state of the art simulations of cosmic structure for¬ 
mation (see, e.g. Agertz et al. 2011; Scannapieco et al. 2012; 
Stinson et al. 2013; Hopkins et al. 2014; Vogelsberger et al. 
2014; Schaye et al. 2015; Khandai et al. 2015). As shown in 
the results overview of Fig. 9, there is no significant differ¬ 
ence between the results obtained for both implementations. 

A possible interpretation of the lack of differences for 
the cosmological runs is that the properties of galaxies and 
their dynamics are already captured sufficiently accurately 


by the standard Arepo implementation. The accuracy im¬ 
provements brought about by the new formulation proposed 
here are either irrelevant for this problem, or are completely 
dominated by first order errors introduced by the sub-grid 
models for radiative cooling, star formation, and feedback. 
The latter is particularly likely as it is now well understood 
that the outcome of galaxy formation simulations depends 
very sensitively on the treatment of highly non-linear feed¬ 
back processes. Thus, in regions of the simulation that are 
crucially shaped by feedback we do not expect our improve¬ 
ments to the code to lead to significant changes of the results. 
Changing this to improve the solution requires implement¬ 
ing and coupling the source terms such that the combined 
system is second order accurate. This may be different in 
regions where hydrodynamics and gravity are the only rele¬ 
vant phenomena. For example, conceivably some differences 
may be found in the approximately hydrostatic atmosphere 
of rich galaxy clusters, where our new formulation may re¬ 
sult in a marginally better representation of turbulence. 


7 CONCLUSIONS 

In this paper, we have discussed in detail two simple modifi¬ 
cations of the cosmological moving-mesh code Arepo, which 
are nevertheless quite important to recover full second-order 
convergence in the L 1 norm for general smooth problems 
with non-trivial mesh motions. One of these changes con¬ 
cerns the time integration, where a Runge-Kutta time inte¬ 
gration scheme is adopted instead of the MUSCL-Hancock 
approach in order to account for changes of the mesh ge¬ 
ometry during a timestep at second order. This does not 
increase the number of mesh-constructions needed, but does 
double the number of required flux computations. The other 
change is the adoption of a more general (and slightly more 
expensive) gradient estimate that retains the necessary ac¬ 
curacy even for large offsets between the mesh-generating 
points and the centers of masses of cells. 

These improvements are most relevant for smooth, pure 
hydrodynamics problems where the influence of self-gravity 
and of complicated source terms is limited. Among the as- 
trophysical problems that fall into this regime and that have 
already been tackled with the original version of Arepo are 
proto-planetary discs (Munoz et al. 2014), cold gas in galac¬ 
tic discs (Smith et al. 2014), and dynamical stellar mergers 
(Pakmor et al. 2013). Such simulations and similar problems 
will benefit in the future from the added accuracy facilitated 
by the improvements proposed here. It is also prudent to 
carry out additional tests to see whether previous simula¬ 
tions were sometimes degradated in a noticeable way by an 
accuracy loss of the code, for example by an unnecessarily 
large error in the conservation of angular momentum. We 
expect such noticeable errors only for problems that evolve 
a system for many dynamical timescales. 

Cosmological simulations of galaxy formation seem to 
be unaffected by the improvements of the hydrodynamical 
moving-mesh scheme proposed here. This is of course re¬ 
assuring as it means that previous results from galaxy for¬ 
mation studies carried out with Arepo, both with zoom-in 
techniques (Marinacci et al. 2013; Pakmor et al. 2014) and 
in large cosmological boxes (Vogelsberger et al. 2014; Genel 
et al. 2014), have not suffered from the convergence rate is- 
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Figure 9. Comparison of two cosmological zoom runs of a galaxy similar to the Milky-Way (as in Marinacci et al. 2013) using either 
the combination of MUSCL-Hancock time integration and improved Green-Gauss gradient estimate, or the combination of Runge-Kutta 
time integration and least squares gradient estimate, respectively. Initial conditions, code configuration and parameters including all 
sub-grid physics are identical for both runs and chosen in line with Marinacci et al. (2013) and Pakmor et al. (2014). The top row shows, 
from left to right, stellar projections of the two runs at z = 0 and the evolution of the star formation rate in the disc of the main galaxy. 
The bottom row shows circularities of stars in the disc, evolution of the total mass of the halo compared to its stellar mass, and the 
average root mean square magnetic field strength in the disc for the two runs. 


sues discussed above. Nevertheless, it is clearly desirable to 
use the improved scheme in the future in this regime as well, 
especially since it offers higher accuracy at an insignificant 
increase of the computational cost. 
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APPENDIX A: YEE VORTEX 
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We choose the parameters as Ti n f = 1, 7 = 1.4, and /3 = 5. 

The setup of a Cartesian mesh is trivial. To organise 
the mesh-generating points on rings around the center of the 
vortex in a regular way, we first compute the width of a ring 
as d r ing = 10/A, where 10 is the size of our domain in one 
dimension and N the linear resolution. We then add mesh¬ 
generating points on circles with radii that are multiples of 
the width of a ring, r r j ng = i x d r i ng . On every circle we 
equidistantly add 2-7T r r in g /d r i ng points. We only add points 
that lie in the computational domain and repeat this until 
none of the points of a new ring are in the domain anymore. 


APPENDIX B: KEPLERIAN DISC 


The setup for the Keplerian disc is similar to the one 
used in Hopkins (2014). We use a computational domain 
of [—2.5, 2.5] x [—2.5, 2.5]. The disc is cold, i.e. the pressure 
support is negligible compared to gravitational force and or¬ 
bital velocity. The initial mesh is again organised on rings 
around the center of the disk. We set the initial surface den¬ 
sity p , velocity v, and specific internal energy u of the cells 
to the values of the continuous fields at the centres of mass 
of the cells. The continuous fields at position r = (*, y) and 
r = I?'| are given as follows. For r < 0.5 and r > 2 we set 


p = 10 

v = 0, 

57 

2 p 

And for 0.5 < r < 2 we adopt 


= —1 x 10 " 



We choose an adiabatic index of 7 = 5/3. The constant 
external gravitational acceleration is given by 


r (r 2 + e 2 ) ’ 


(Bl) 


For the setup of the Yee vortex we follow Yee et al. (2000). 
The mesh is defined on the domain [—5, 5] x [—5, 5]. In the 
initial conditions, we set density p, velocity v, and specific 
internal energy u of the cells to the value of the continuous 
fields at the centres of mass of the cells. The continuous 
fields at position r = (x, y) are given by 


where e = 0.25 for r < 0.25 and t = 0 everywhere else. 

This paper has been typeset from a J’pN/bHjrX file prepared by 
the author. 
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